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Abstract 

In this short note, we present a new technique to accelerate the convergence of a FFT- 
based solver for numerical homogenization of complex periodic media proposed by 
Moulinec and Suquet [1]. The approach proceeds from discretization of the governing 
integral equation by the trigonometric collocation method due to Vainikko [2], to 
give a linear system which can be efficiently solved by conjugate gradient methods. 
Computational experiments confirm robustness of the algorithm with respect to its 
internal parameters and demonstrate significant increase of the convergence rate for 
problems with high-contrast coefficients at a low overhead per iteration. 

Key words: numerical homogenization; FFT-based solvers; trigonometric 
collocation method; conjugate gradient solvers 
PACS: 72.80.Tm; 72.80.-m 



1. Introduction 

A majority of computational homogenization algorithms rely on the solution of 
the unit cell problem, which concerns the determination of local fields in a repre- 
sentative sample of a heterogeneous material under periodic boundary conditions. 
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Currently, the most efficient numerical solvers of this problem are based on dis- 
cretization of integral equations. In the case of particulate composites with smooth 
bounded inclusions embedded in a matrix phase, the problem can be reduced to 
internal interfaces and solved with remarkable accuracy and efficiency by the fast 
multipole method, see [3, and references therein]. An alternative method has been 
proposed by Suquet and Moulinec [1] to treat problems with general microstructures 
supplied in the form of digital images. The algorithm is based on the Neumann series 
expansion of the inverse to an operator arising in the associated Lippmann-Schwinger 
equation and exploits the Fast Fourier Transform (FFT) to evaluate the action of 
the operator efficiently. 

The major disadvantage of the FFT-based method consists in its poor conver- 
gence for composites exhibiting large jumps in material coefficients. To overcome 
this difficulty, Eyre and Milton proposed in [4] an accelerated scheme derived from 
a modified integral equation treated by means of the series expansion approach. In 
addition, Michel et al. [5] introduced an equivalent saddle-point formulation solved 
by the Augmented Lagrangian method. As clearly demonstrated in a numerical 
study by Moulinec and Suquet [6], both methods converge considerably faster than 
the original variant; the number of iterations is proportional to the square root of 
the phase contrast instead of the linear increase for the basic scheme. However, this 
comes at the expense of increased computational cost per iteration and the sensitivity 
of the Augmented Lagrangian algorithm to the setting of its internal parameters. 

In this short note, we introduce yet another approach to improve the conver- 
gence of the original FFT-based scheme [1] based on the trigonometric colloca- 
tion method [7] and its application to the Helmholtz equation as introduced by 
Vainikko [2]. We observe that the discretization results in a system of linear equations 
with a structured dense matrix, for which a matrix-vector product can be computed 
efficiently using FFT, cf. Section 2. It is then natural to treat the resulting system 
by standard iterative solvers, such as the Krylov subspace methods, instead of the 
series expansion technique. In Section 3, the potential of such approach is demon- 
strated by means of a numerical study comparing the performance of the original 
scheme and the conjugate- and biconjugate- gradient methods for two-dimensional 
scalar electrostatics. 

2. Methodology 

In this section, we briefly summarize the essential steps of the trigonometric 
collocation-based solution to the unit cell problem by adapting the original exposition 
by Vainikko [2] to the setting of electrical conduction in periodic composites. In what 
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follows, a, a and A denote scalar, vector and second-order tensor quantities with 
Greek subscripts used when referring to the corresponding components, e.g. A a p. 
Matrices are denoted by a serif font (e.g. A) and a multi-index notation is employed, 
in which R N with N = (N 1 ,...,N d ) represents R-^i*-*^ an d A fc stands for the 
(ki, . . . , kd)-th element of the matrix A G R . 

2.1. Problem setting 

We consider a composite material represented by a periodic unit cell y = Ylt=i(~^a, Y<x) C 
M. d . In the context of linear electrostatics, the associated unit cell problem reads as 

Vxe(a;) = 0, V-j(x) = 0, j(x) = L(x) ■ e(x), x G y (1) 

where e is a ^-periodic vectorial electric field, j denotes the corresponding vector of 
electric current and L is a second-order positive-definite tensor of electric conductiv- 
ity. In addition, the field e is subject to a constraint 

e ° = W\J,« x)ix < (2) 

where e° denotes a prescribed macroscopic electric field and |3^| represents the d- 
dimensional measure of y. 

Next, we introduce a homogeneous reference medium with constant conductivity 
L , leading to a decomposition of the electric current field in the form 

j(x) = L°-e(x) + 5L(x)-e(x), 5L(x) = L(x) - L°. (3) 

The original problem (l)-(2) is then equivalent to the periodic Lippmann-Schwinger 
integral equation, formally written as 

e(x) + J r°(x - y) ■ (sL(y) ■ e{yj) dy = e°, x G y, (4) 

where the T operator is derived from the Green's function of the problem (l)-(2) 
with L(x) = L° and e° = 0. Making use of the convolution theorem, Eq. (4) attains 
a local form in the Fourier space: 

e(k) = l n Ji |le °' fc = °' (5) 
\ -P {k)-{5L-e){k), keZ d \{0}, 

where f(k) denotes the Fourier coefficient of f(x) for the fc-th frequency given by 

f(k) = [ f(x)<p- k (x)dx, <p k (x) = |^|"5exp ( m V ) , (6) 
J y \ a=i Ya J 
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"i" is the imaginary unit and 

f 0, fe = 0, 

Here, we refer to [4, 8] for additional details. 

2.2. Discretization via trigonometric collocation 

Numerical solution of the Lippmann-Schwinger equation is based on a discretiza- 
tion of a unit cell y into a regular periodic grid with N±X ■ ■ ■ X N<i nodal points and 
grid spacings h = (2Yi/iVi, . . . , 2Yd/Nd). The searched field e in (4) is approximated 
by a trigonometric polynomial e N in the form (cf. [2]) 

e(x) « e N (x) = e(k)v>k{*), (8) 
kez N 

where N = (Ni, . . . , Nj), e designates the Fourier coefficients defined in (6) and 

z N = {fcez d :4<t<y' a = 1 ''"4 (9) 

We recall, e.g. from [2], that the a-th component of the trigonometric polynomial 
expansion admits two equivalent finite-dimensional representations. The first one 
is based on a matrix e Q G C"^ of the Fourier coefficients of the a-th component and 
equation (8) with e a (k) = e k . Second, the data can be entirely determined by 
interpolation of nodal values 



e?(x) 



X)eS^(«), a = l,...,d (10) 
kez N 



where e a G M. N is a matrix storing electric field values at grid points, = (x k ) 
is the corresponding value at the k-th node with coordinates x k = (kih\, . . . , kjid) 
and basis functions 



) = \N\-^e W L± m J^- 2 M\ (11) 
me z« I o=i Vla iVo/ J 

satisfy the Dirac delta property (p™ (x m ) = Smk with |JV| = Ylt=i N a - Both repre- 
sentations can be directly related to each other by 

e Q Fe Q , 6 a F e Q , (12) 
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where the Vandermonde matrices F 6 £^ NxN an( j F 1 6 £^ NxN 

pfem 



(F- 1 )^ = l^jiY^exp ]>>i^ 



v a=l 



(13) 
(14) 



implement the forward and inverse Fourier transform, respectively, e.g. [9, Sec- 
tion 4.6]. 

The trigonometric collocation method is based on the projection of the Lippmann- 
Schwinger equation (4) to the space of the trigonometric polynomials of the form 
{^2k& N c feVfe> c k w, cf. [7, 2]. In view of Eq. (10), this is equivalent to the 
collocation at grid points, with the action of T° operator evaluated from the Fourier 
space expression (5) converted to the nodal representation by (12)2. The resulting 
system of collocation equations reads 



:i + B)e = e c 



(15) 



where e 6 R dxN and e° 6 



pdxN 



store the corresponding solution and of the macro- 



scopic field, respectively. Furthermore, I is the d x d x N x N unit matrix and the 
non-symmetric matrix B can be expressed, for the two-dimensional setting, in the 
partitioned format as 



B 



F 1 






F- 1 



r° 
1 11 

1 21 



r° 
1 12 

r° 

1 22 



F 
F 



5L U 5l 12 

51-21 SL 2 2 



(16) 



with an obvious generalization to an arbitrary dimension. Here, V 



it holds 





a/3 



iNxN 



and 



are diagonal matrices storing the corresponding grid values, for which 



fefe 



a/3 



X 



a, (3 



,dzndkeZ N . (17) 



2.3. Iterative solution of collocation equations 

It follows from Eq. (16) that the cost of the multiplication by B or by B T is 
driven by the forward and inverse Fourier transforms, which can be performed in 
0(|iV| log I A?" |) operations by FFT techniques. This makes the resulting system (15) 
ideally suited for iterative solvers. 
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In particular, the original Fast Fourier Transform-based Homogenization (FFTH) 
scheme formulated by Moulinec and Suquet in [1] is based on the Neumann expansion 
of the matrix inverse (I + B) -1 , so as to yield the m-th iterate in the form 

rn 

e (m) =J2(-B) j e°. (18) 

3=0 

Convergence of the series (18) was comprehensively studied in [4, 8], where it was 
shown that the optimal rate of convergence is achieved for 

L° = Amin + Anxax I, (19) 

with A m i n and A max denoting the minimum and maximum eigenvalues of L{x) on y 
and I being the identity tensor. 

Here, we propose to solve the non-symmetric system (15) by well-established 
Krylov subspace methods, in particular, exploiting the classical Conjugate Gradi- 
ent (CG) method [10] and the biconjugate gradient (BiCG) algorithm [11]. Even 
though that CG algorithm is generally applicable to symmetric and positive-definite 
systems only, its convergence in the one- dimensional setting has been proven by 
Vondfejc [12, Section 6.2]. A successful application of CG method to a generalized 
Eshelby inhomogeneity problem has also been recently reported by Novak [13] and 
Kanaun [14]. 

3. Results 

To assess the performance of the conjugate gradient algorithms, we consider a 
model problem of the transverse electric conduction in a square array of identical cir- 
cular particles with 50% volume fraction. A uniform macroscopic field e° = (1, 0) is 
imposed on the corresponding single-particle unit cell, discretized by N = (255, 255) 
nodes 1 and the phases are considered to be isotropic with the conductivities set to 
L = I for the matrix phase and to L = gl for the particle. 

The conductivity of the homogeneous reference medium is parameterized as 



•o, 



u) = (1 - w + gu) I, (20) 



1 Note that the odd number of discretization points is used to eliminate artificial high-frequency 
oscillations of the solution in the Fourier space, as reported in [15, Section 2.4]. 
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where u = 0.5 corresponds to the optimal convergence of FFTH algorithm (19). 
All conjugate gradient-related results have been obtained using the implementa- 
tions according to [16] and referred to as Algorithm 6.18 (CG method) and Algo- 
rithm 7.3 (BiCG scheme). Two termination criteria are considered. The first one is 
defined for the m-th iteration as [15] 



and provides the test of the equilibrium condition (1)2 in the Fourier space. An 
alternative expression, related to the standard residual norm for iterative solvers, 
has been proposed by Vinogradov and Milton in [8] and admits the form 



with the additional L° term ensuring the proportionality to (21) at convergence. 
From the numerical point of view, the latter criterion is more efficient than the 
equilibrium variant, which requires additional operations per iteration. From the 
theoretical point of view, its usage is justified only when supported by a convergence 
result for the iterative algorithm. In the opposite case, the equilibrium norm appears 
to be more appropriate, in order to avoid spurious non-physical solutions. 

3.1. Choice of reference medium and norm 

Since no results for the optimal choice of the reference medium are known for 
(Bi)CG-based solvers, we first estimate their sensitivity to this aspect numerically. 
The results appear in Fig. 1(a), plotting the relative number of iterations for CG and 
BiCG solvers against the conductivity of the reference medium parameterized by u>, 
recall Eq. (20). 

As expected, both CG and BiCG solvers achieve a significant improvement over 
FFTH method in terms of the number of iterations, ranging from 50% for a mildly- 
contrasted composite down to 2% for g = 10 4 . Moreover, contrary to all other 
available methods, the number of iterations is almost independent of the choice of 
the reference medium. We also observe, in agreement with results by [12, Section 6.2] 
for the one-dimensional setting, that CG and BiCG algorithms generate identical 
sequences of iterates; the minor differences visible for u > 1 or g = 10 4 can be 
therefore attributed to accumulation of round-off errors. These conclusions hold for 
both equilibrium- and residual-based norms, which appear to be roughly proportional 




(21) 



(m) _ 




< 6 



(22) 
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Figure 1: (a) Relative number of iterations as a function of the reference medium parameter ui and 
(b) ratio between residual- and equilibrium-based norms at convergence for r\ x termination condition 
with tolerance e = 10~ 4 . 



for the considered range of the phase contrasts, cf. Fig. 1(b). Therefore, the residual 
criterion (22) will mostly be used in what follows. 

In Fig. 2, we supplement the comparison by considering the total CPU time 
required to achieve a convergence. The data indicate that the cost of one iteration 
is governed by the matrix- vector multiplication, recall Eq. (16): the overhead of CG 
scheme is about 10% with respect to FFTH method, while the application of BiCG 
algorithm, which involves B and B T products per iteration [11], is about twice as 
demanding. As a result, CG algorithm significantly reduces the overall computational 
time in the whole range of contrasts, whereas a similar effect has been reported for 
the candidate schemes only for g > 10 3 , cf. [6]. 

3.2. Influence of phase contrast 

As confirmed by all previous works, the phase contrast g is the critical parameter 
influencing the convergence of FFT-based iterative solvers. In Fig. 3, we compare the 
scaling of the total number of iterations with respect to phase contrast for CG and 
FFTH methods, respectively. The results clearly show that the number of iterations 
grows as ^/g instead of the linear increase for FFTH method. This follows from error 
bounds 

4 m) <l m vi 0) , 7 FFTH = ^, 7 CG = ^|. (23) 

0+1 a/0+ 1 
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Figure 2: Relative CPU time t of CG and BiCG solvers plotted against the conductivity parameter 
oj for ?7 r -based termination condition with tolerance e = 10~ 4 . 

The first estimate was proven in [4], whereas the second expression is a direct conse- 
quence of the condition number of matrix B being proportional to q and a well-known 
result for the conjugate gradient method, e.g. [16, Section 6.11.3]. The CG-based 
method, however, failed to converge for the infinite contrast limit. Such behavior is 
equivalent to the Eyre-Milton scheme [4]. It is, however, inferior to the Augmented 
Lagrangian algorithm, for which the convergence rate improves with increasing p and 
the method converges even as p — > oo. Nonetheless, such results are obtained for 
optimal, but not always straightforward, choice of the parameters [5]. 

3.3. Convergence progress 

The final illustration of the CG-based algorithm is provided by Fig. 4, displaying 
a detailed convergence behavior for both low- and high-contrast cases. The results 
in Fig. 4(a) correspond well with estimates (23) for both residual and equilibrium- 
based norms. Influence of a higher phase contrast is visible from Fig. 4(b), plotted 
in the full logarithmic scale. For FFTH algorithm, two regimes can be clearly dis- 
tinguished. In the first few iterations, the residual error rapidly decreases, but the 
iterates tend to deviate from equilibrium. Then, both residuals are simultaneously 
reduced. For CG scheme, the increase of the equilibrium residual appears only in 
the first iteration and then the method rapidly converges to the correct solution. 
However, its convergence curve is irregular and the algorithm repeatedly stagnates 
in two consecutive iterations. Further analysis of this phenomenon remains a subject 
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phase contrast p 



Figure 3: Total number of iterations n plotted against phase contrast g; the reference medium 
corresponds to for to = 0.5 and tolerance s is related to rj r norm. 



of future work. 
4. Conclusions 

In this short note, we have presented a conjugate gradient-based acceleration of 
the FFT-based homogenization solver originally proposed by Moulinec and Suquet [1] 
and illustrated its performance on a problem of electric conduction in a periodic 
two-phase composite with isotropic phases. On the basis of obtained results, we 
conjecture that: 

i) the non-symmetric system of linear equations (15), arising from discretization 
by the trigonometric collocation method [2], can be solved using the standard 
conjugate gradient algorithm, 

ii) the convergence rate of the method is proportional to the square root of the 
phase contrast, 

iii) the methods fails to converge in the infinite contrast limit, 

iv) contrary to available improvements of the original FFT-solver [4, 5], the cost 
of one iteration remains comparable to the basic scheme and the method is 
insensitive to the choice of auxiliary reference medium. 
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(a) (b) 

Figure 4: Convergence progress of CG and FFTH methods for (a) g = 10 1 and (b) g = 10 3 as 
quantified by rj e and rj r norms; reference medium corresponds to u — 0.5 and the dot-and-dahsed 
curves indicate the convergence rates (23). 

The presented computational experiments provide the first step towards further 
improvements of the method, including a rigorous analysis of its convergence prop- 
erties, acceleration by multi-grid solvers and preconditioning and the extension to 
non-linear problems. 
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